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We compute the QGP suppression of T(ls), T(2s), T(3s), Xbi, and Xb 2 states in 

,ysNN = 2.76 TeV Pb-Pb collisions. Using the suppression of each of these states, we estimate 
the inclusive Raa for the T(ls) and T(2s) states as a function of Vpart, y, and pt including the 
effect of excited state feed down. We find that our model provides a reasonable description of 
preliminary CMS results for the Vpart-, y-, and pT-dependence of Raa for both the T(ls) and 
T(2s). Comparing to our previous model predictions, we find a flatter rapidity dependence, thereby 
reducing some of the tension between our model and ALICE forward-rapidity results for T(ls) 
suppression. 
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The relativistic heavy-ion collision experiments being 
carried out at Brookhaven National Laboratory’s Rela¬ 
tivistic Heavy Ion Collider (RHIC) and CERN’s Large 
Hadron Collider (LHC) study the behavior of matter at 
extreme temperatures and densities. The goal of these 
experiments is to generate a deconfined state of nuclear 
matter called a quark-gluon plasma (QGP) and to study 
its properties in detail. Based on hydrodynamic fits 
to particle production, LHC ^/SNN = 2.76 TeV colli¬ 
sions generate QGP initial temperatures on the order of 
To ~ 500 — 600 MeV [1, 2]. At such high temperatures 
light hadronic states are disassociated and the equation 
of state of nuclear matter is well-described by a gas of 
quark and gluon quasiparticles [3, 4]. In the transition 
region between hadronic matter and a proper QGP, the 
system is composed of liberated quarks and gluons plus 
a small admixture of heavy bound states. Although light 
hadronic states disassociate around the pseudo-critical 
temperature for the quark-hadron transition, Tc ~ 165 
MeV, bottomonia, for example, may survive up to tem¬ 
peratures on the order of T ~ 600 MeV ~ 4Tc [5]. Due 
to mass/binding-energy ordering of the quarkonium spec¬ 
trum, one expects that there will be an approximate se¬ 
quential disassociation, with lighter states “melting” be¬ 
fore heavier states and excited states melting before their 
respective ground states [6]. 

In this paper, we focus on the suppression of bottomo¬ 
nia in ytsNN = 2.76 TeV Pb-Pb collisions. The bene¬ 
fits of working with heavy quarks are that heavy quark 
bound states are dominated by short distance physics, 
their binding energies are much smaller than the quark 
mass rriQ Aqcd (Q = c, 6), and their sizes are much 
larger than IjmQ. As a result, they can be treated using 
effective field theory methods. In the heavy quark limit, 
one finds that a potential-based non-relativistic effective 
field theory, pNRQGD, can be used to calculate the mass 
spectrum, decay rates, etc. of heavy quark bound states 
[7-11]. In addition, pNRQCD allows for the systematic 


inclusion of relativistic corrections. Using pNRQCD po¬ 
tential models, the vacuum spectrum of all bottomonium 
states can be reproduced to within less than one percent 
using a Cornell potential plus spin-spin and spin-orbit 
interactions [12, 13]. 

The use of potential models to describe quarkonium 
suppression has a long history, starting with the seminal 
works of Karsch, Matsui, Mehr, and Satz [14, 15] who 
predicted that quarkonium production would be sup¬ 
pressed in heavy-ion collisions due to Debye-screening in 
a deconfined QGP. Using such non-relativistic potential 
models, there have been studies of quarkonium spectral 
functions and mesonic current correlators, see e.g. [16- 
23] . There have also been lattice QCD calculations of the 
quarkonium spectral function [24-33]. Compared to the 
standard Debye-screened potential models used in early 
calculations, systematic analysis of the heavy quark po¬ 
tential in the QGP showed that the potential is complex¬ 
valued, with the imaginary part of a state’s energy being 
related to the thermal width of the state [34]. 

In the bottom sector, potential model calculations in¬ 
dicate that the T(ls), T(2s), and T(3s) can survive up 
to temperatures T ~ 593, 228, 172 MeV, respectively 
[5]. At these temperatures, the in-medium width of the 
state becomes on the order of the real part of its binding 
energy, and the bound state quickly disappears from the 
spectrum. However, even below this disassociation point, 
quarkonia also decay due to in-medium interactions. For 
the T(ls), the in-medium width approaches 100 MeV at 
3Tc [35]. At this temperature, the T(ls) in-medium half- 
life is on the order 2 fm/c. Since this is also the timescale 
over which the QGP evolves hydrodynamically, one needs 
accurate and reliable modeling of the background evolu¬ 
tion in order to make reliable predictions for quarkonium 
suppression in heavy-ion collisions. 

In this paper, we provide an update to the model used 
in Refs. [36, 37] to: (1) extend the background evolu¬ 
tion to full (3-|-l)D anisotropic hydrodynamics (aHydro) 
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with a rapidity profile consistent with experimentally- 
observed particle multiplicity distributions; (2) update 
the mixing fractions to recent updated values determined 
via fits to ATLAS, CMS, and LHCb results for T and Xb 
production in p-p collisions [38]; (3) correct the proba¬ 
bility weight-function used for centrality averaging in or¬ 
der to match the experimental procedure. We compare 
the updated model predictions with recently reported re¬ 
sults on T suppression in Pb-Pb collisions from both the 
CMS [39] and ALICE [40] collaborations. We find that, 
with the improvements listed above, the original model of 
Refs. [36, 37] gives a reasonable description of the A^part-, 
?/-, and pT-dependence of T(ls) and T(2s) suppression. 

Methodology. For both the potential and dynamical 
equations used below, we assume that the effective lo¬ 
cal rest frame (LRF) one-particle distribution function 
for the particles comprising the QGP is of the form 

f{p,x) = feq(^\/pT + [I + ^{x)\pI ! K{x'^ , ( 1 ) 


where — 1 < ^(a;) < cx) is the local spheroidal momentum- 
space anisotropy parameter and A(a;) is the local trans¬ 
verse temperature. This form takes into account the 
difference between the transverse and longitudinal pres¬ 
sures, which is the most important viscous correction 
generated in heavy-ion collisions. 

As mentioned above, it is now understood that the 
heavy quark potential in the QGP has both real and 
imaginary parts, V = 3fi[F] -|-f5[P]. We use the internal- 
energy-based potential specified originally in Ref. [37]. 
In the model, the real part of the potential is obtained 
from the internal energy of the heavy quark/anti-quark 
system [56]. The resulting real part of the potential is 
given by [37] 
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where mb = 4.7 GeV, a = 0.385, ct = 0.223 GeV^ [41], 
and the last term is a temperature- and spin-independent 
finite-quark-mass correction taken from Ref. [42] . In this 
expression, fi = Gi^,d)mD [37, 43, 44] is the anisotropic 
Debye mass, where ^ is a function which depends on 
the degree of plasma momentum-space anisotropy the 
angle of the line connecting the quark-antiquark pair 
with respect to the beamline direction 0, and mu = 
1.4-y/l -I- Nf/6 gsT is the isotropic leading-order Debye 
mass adjusted by a factor of 1.4 in order to take into 
account higher-order corrections determined via lattice 
calculations [45]. Note that, in the limit ^ 0, one has 

G = 1 and the real part of the potential above reduces 
to the internal energy derived from the original Karsch- 
Mehr-Satz potential [15]. 


The imaginary part of the potential 3[F] is obtained 
from a leading-order perturbative calculation performed 
in the small-^ limit [34, 46, 47] 

3[F] = -OsCpT 

-CbPiir/mD,0) , (3) 

where (p, '01 j 02 are special functions which can be 

expressed in terms of the Meijer G-function. We solve 
the 3D Schrodinger equation with the potential above to 
obtain the real and imaginary parts of the binding energy 
as a function of ^ and A [35] . The imaginary part of the 
binding energy is then used to obtain the width of each 
state 


r(T, x_L,0 


23[£;bi„d] 5i[F;bi„d] > 0 

7dis ^[Ebind]<0 ’ 


with 7dis being the effective decay rate for unbound 
states, which we take to be 10 GeV [57]. It is implicitly 
understood that Fbind, and hence T, are local quantities 
that depend on r = xj_, and = tanh“^(z/t) 

through the (3-|-l)D evolution of the transverse temper¬ 
ature A, local momentum-space anisotropy .f, and asso¬ 
ciated flow velocities. For this purpose, we use (3-|-l)D 
anisotropic hydrodynamics (aHydro) [48-50]. 

The (3-|-l)D aHydro code used provides the spatiotem- 
poral evolution of ^ and A. The widths obtained from 
solution of the 3D Schrodinger equation are then inte¬ 
grated and exponentiated to compute the relative num¬ 
ber of states remaining at a given proper time. Integrat¬ 
ing the instantaneous local decay rate T over proper-time, 
one obtains 


C = 0(t/ - Tform) / dr r(T, X^, ■? = 1/) , (5) 

dmax(rform,To) 

where b is the impact parameter, Tform = TRorm = 
'^TTfoj.^/M where M is the mass of the state, and is 
the formation time of the state at rest. For the rest frame 
formation times, we assume that they are roughly pro¬ 
portional to the inverse vacuum binding energy of each 
of the states [51]. For the T(ls), T(2s), T(3s), x^i, and 
Xb 2 states, we use = 0.2, 0.4, 0.6, 0.4, and 0.6 fm/c, 
respectively [58]. 

We take the initial proper time tq for hydrody¬ 
namic evolution to be tq = 0.3 fm/c and the ini¬ 
tial central temperature for central collisions to be 
Tq € {552,546,544} MeV for shear viscosity to entropy 
density ratios 47rrj/s € {1,2,3}, with the values tuned 
in order to keep the final charged particle multiplicity 
fixed. The final time ry appearing in Eq. (5) is self- 
consistently determined from the aHydro simulation as 
the proper time when local effective temperature be¬ 
comes less than the transition temperature. At this ef¬ 
fective temperature, plasma screening effects are assumed 
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to decrease rapidly due to the transition to the hadronic 
phase with the widths of the states becoming approxi¬ 
mately equal to their vacuum widths [59]. For the aHy- 
dro initial conditions, we use a smooth linear combina¬ 
tion (Kbinary = 0.145) of Glauber wounded-nucleon and 
binary collision scaling to set the initial energy density 
profile in the transverse plane. The inelastic cross-section 
is taken to be ctjvat = 62 mb. In the spatial rapidity di¬ 
rection, we use a boost-invariant plateau at central ra¬ 
pidities with Gaussian-tails consistent with limited frag¬ 
mentation at large rapidity [52] 


/(<7) = exp 
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( 6 ) 


with Ac = 2.5 and = 1.4 fitted to reproduce the exper¬ 
imental pseudorapidity distribution of charged particles. 

In order to compare to the experimental results, we 
then (a) perform a weighted average over the transverse 
plane and (b) implement any cuts on centrality, pT, and 
rapidity necessary. For the spatial average, the proba¬ 
bility distribution function for bottomonium production 
is taken to be proportional to the local number den¬ 
sity of plasma partons n(xj_,c), i.e. RAAiPr^y^b) = 
(/xi n(x_L,c). For imple¬ 

menting pt cuts, we assume that the pt probability dis¬ 
tribution function is proportional to E^'^. For imple¬ 
menting cuts in rapidity, we use a flat distribution func¬ 
tion. After implementing the appropriate cuts on px and 
y, we obtain i?AA(fc)- We then convert b to centrality C 
using the Glauber formalism and integrate over the ap¬ 
propriate centrality cuts using a probability distribution 
function proportional to e“W20^ where 0 < C < 100. 
This probability distribution function takes into account 
the increased particle production that occurs in central 
collisions and its form is taken from fits to experimentally 
observed centrality distributions [53]. 

The procedure outlined above gives the “raw” suppres¬ 
sion factors for each state. In order to account for post- 
QGP feed down of excited states for the T(ls), we use 
pr-averaged feed down fractions obtained recently from a 
compilation of p-p data available from ATLAS, GMS, and 
LHCb which gives //* = {0.618,0.105,0.02,0.207,0.05} 
for the T(ls), T(2s), T(3s), Xbi, and Xb 2 to T(ls) feed 
down fractions, respectively [38]. For the T(2s), we as¬ 
sume that /f® = (0.5, 0.5} for the T(2s),T(3s) to T(2s) 
feed down fractions, respectively [54]. 

Results. In Fig. 1, we show the raw Raa for the five 
states considered as a function of A^part for the case that 
Airyjs = 1. As can be seen from this hgure, there is a 
sequential suppression of the states, however, there are 
no thresholds visible as originally predicted by sequen¬ 
tial suppression [6]. The lack of thresholds is due to (1) 
averaging over the full temperature distribution in the 
transverse plane where the QGP is hotter in the center 
and colder as one moves towards the edges and (2) the 



FIG. 1: (Color online) Raw Raa as a function of Apart. For 
this figure we assumed ATvyjs = 1. These curves do not in¬ 
clude the effect of excited state feed down. 



FIG. 2: (Color online) Inclusive Raa for the T(ls) and T(2s) 
as a function of Apart. 


continuous decays of the various states prior to their dis- 
association point. Although, we only show results as a 
function of Apart, the model provides the full Apart-, Pt-, 
and y-dependence of Raa for each of the states. 

By constructing a linear combination of the raw Raa 
for each state, we obtain the inclusive Raa for the states. 
The result of performing this procedure for the T(ls) and 
T(2s) is plotted in Figs. 2-4. In these three hgures, each 
set of three lines corresponds to Airyjs G {1,2,3}. In 
Fig. 2, we compare our results to recently reported pre¬ 
liminary data from the GMS collaboration [39]. As can 
be seen from this figure, our model does a good job for 
both the T(ls) and T(2s) states. There is, however, some 
tension with the lowest Apart point for Based on 

the comparison of the model predictions with CMS pre¬ 
liminary data for rJ}a^\ the data seem to prefer small 
shear viscosities in the range 1 < Ekt]! s < 2. The R^^pf^ 
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FIG. 3: (Color online) Inclusive Raa for the T(ls) and T(2s) 
as a function of y. 



FIG. 4: (Color online) Inclusive Raa for the T(ls) and T(2s) 
as a function of pr- 


data does not seem to provide a tight constraint on 77 /s 
at this point in time. 

In Fig. 3, we show our results as a function of rapidity 
and, once again, we compare with the new CMS prelimi¬ 
nary data. We also include the result obtained by 

the ALICE collaboration at forward rapidities as open 
circles [40]. Although our model does a reasonable job 
in reproducing the trends seen in the CMS preliminary 
data, there is still some lingering tension with the AL¬ 
ICE forward results. We note, however, that compared 
with earlier predictions made in Ref. [54], our model re¬ 
sults are now much closer to the ALICE data. This is 
due solely to the change in the way we perform the cen¬ 
trality averaging. In the past, we used a flat probability 
distribution as a function of centrality, which does not 
conform to the procedure used to compute the centrality- 
averaged results by the experiments, where they simply 
average over the particles detected in each centrality bin. 
With the updated probability distribution function, the 


centrality-averaged results are much closer to those ob¬ 
tained in central collisions. 

Finally, in Fig. 4 we show our results as a function of 
Pt compared to CMS preliminary data. The flatness of 
Raa as a function of pr was a prediction contained in 
the original model [36, 37] and is due to the fact that, in 
the model, the bottomonia spectra are assumed to be un¬ 
affected due to the lack of thermalization of these states 
because of their large masses. The slow increase in Raa 
as a function of pr stems solely from the effect of time- 
dilation of the formation times of the states. Comparing 
to the CMS preliminary results for aJa 
data seem to, once again, prefer small values of 77 /s. For 
the rJa model seems to under predict the amount 

of suppression seen in the CMS preliminary data, how¬ 
ever the overall magnitude and weak dependence on pj- 
predicted by the model seems to be in reasonable agree¬ 
ment with the data. 

Conclusions. In this paper we presented an update to 
our model predictions for the QGP-induced suppression 
of bottomonia states at LHC energies [36, 37]. The po¬ 
tential model itself is exactly the same as used in previ¬ 
ously published results, however, we have (1) upgraded 
the aHydro code to (3-|-l)D in order to have a more real¬ 
istic model of the background evolution (2) updated the 
mixing fractions determined from recent ATLAS, CMS, 
and LHCb measurements, and (3) corrected our method 
for performing centrality averaging. 

As can be seen from the results presented herein, the 
original internal-energy-based model of Refs. [36, 37] 
seems to do a reasonable job describing the A^part-, 2/-, 
andpT-dependence of CMS preliminary results for rJ}a^^ 
and At forward rapidities, there is still some ten¬ 

sion with the ALICE rJa data, however, with the fix 
to the centrality-averaging procedure, the discrepancy is 
no longer as dramatic. Because of this, there is now some 
hope that the additional suppression at forward rapidi¬ 
ties could be explained by cold-nuclear matter effects. 

On the positive side, it seems that for central rapidi¬ 
ties (77 < 2) the data are consistent with bottomonia 
suppression due to the creation of a deconfined QGP 
with a shear viscosity to entropy density ratio roughly 
between l/(47r) and 2/(47r). These values are consistent 
with those obtained via analysis of the collective flow 
coefficients, thereby providing further evidence that the 
QGP created in relativistic heavy ion collisions behaves 
like a nearly perfect fluid. 
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